THE SDA METHOD FOR NUMERICAL SOLUTION OF LUR'E 

EQUATIONS 

FEDERICO POLONF AND TIMO REISt 

Abstract. Wc introduce a numerical method for the numerical solution of the so-called Lur'e 
matrix equations that arise in balancing-related model reduction and linear-quadratic infinite time 
horizon optimal control. Based on the fact that the set of solutions can be characterized in terms of 
deflating subspaces of even matrix pencils, an iterative scheme is derived that converges linearly to 
the maximal solution. 
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1. Introduction. For given matrices A,Q ^ i£n,n .^^ith Q ~ Q* and B,C € 

C"'™, R e c™'™, we consider Lur'e equations 

A*X + XA + Q = K*K, 

XB + C = K*L, (1.1) 
R = L*L, 

that have to be solved for {X, K, L) e C"-" x C*'-" x C^-™ with X = X* and p as 
small as possible. Equations of type (|l.ip were first introduced by A.I. Lur'e [HT] 
in 1951 (see |S] for an historical overview) and play a fundamental role in systems 
theory, e.g. since properties like dissipativity of linear systems can be characterized via 
their solvability [TJ [2l 131 143] . This type of equations moreover appears in the infinite 
time horizon linear-quadratic optimal control problem [T2j [13l [TH |44l |45] , spectral 
factorization [TU [15] as well as in balancing-related model reduction O EH [33l [34l [37] . 
In the case where R is invertible, the matrices K and L can be eliminated by obtaining 
the algebraic Riccati equation 

A*X + XA- {XB + C)R~\XB + 0)* +Q = 0. (1.2) 

Whereas this type is well-explored both from an analytical and numerical point of view 
[Ml [HI [7], the case of singular R has received comparatively little attention. However, 
the singularity of R is often a structural property of the system to be analyzed |36| 
and can therefore not be avoided by arguments of genericity. 

Two approaches exist for the numerical solution of Lur'e equations with (possibly) 
singular R. The works |26l [S] present an iterative technique for the elimination of 
variables corresponding to kcr R. After a finite number of steps this leads to a Riccati 
equation. This also gives an equivalent solvability criterion that is obtained by the 
feasibility of this iteration. The most common approach to the solution of Lur'e 
equations is the slight perturbation of R by e/,„ for some e > 0. Then by using the 
invertibility of R+sI, the corresponding perturbed Lur'e equations are now equivalent 
to the Riccati equation 

A*X^ + X^A - {XB + C){R + eiy\X^B + C)* + Q = 0. (1.3) 
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It is shown in [571 ISH] that certain corresponding solutions X,, then converge to a 
solution of (|l.ip . 

Whereas the first approach has the great disadvantage that it relies on successive 
nullspace computations (which may be arbitrarily an ill-conditioned numerical prob- 
lem), the big problem of the perturbation approach is that there exist no estimates 
for the perturbation error \\X — X^W and, furthermore, the numerical condition of the 
Riccati equation (|1.3p increases drastically as e tends to 0. 

From a theoretical point of view, Lur'e equations have been investigated in |15l 
The solution set is completely characterized in |35| via the consideration of the matrix 
pencil 



This pencil has the special property of being even, that is £ is skew-Hermitian and 
A is Hermitian. Solvability of (jl.ip is characterized via the eigenstructure of this 
pencil. It is furthermore shown that there exists some correspondence to deflating 
subspaces of ()1.4|) . That is, a generalization of the concept of invariant subspaces to 
matrix pencils (TB]. Under some slight additional conditions of the pair {A,B), it is 
shown in |35] that there exists a so-called maximal solution X. In this case, maximal 
means that X is, in terms of semi-definiteness, above all other solutions of the Lur'e 
equations. This solution is of particular interest in optimal control as well as in model 
reduction. 

Based on these theoretical results of [35], we will set up an iterative scheme that 
converges linearly to the maximal solution. 

The paper is organized as follows. Section [2| introduces the notation and contains 
some required control and matrix theoretic background, in particular a normal form 
for even matrix pencils is introduced. In Section|3]we briefly present some results from 
[35] which connect the spectral properties of the even pencil ()1.4|) to the solvability and 
the solutions of the Lur'e equations. An outline of the structured doubling algorithm 
for the solution of control problems is presented in Section (4] Section [5] contains a 
method for transforming a matrix pencil in the form required by the SDA. In Section[6l 
we apply this method to the singular control problem to obtain a reduced discrete-time 
pencil associated to a Lur'e equation. Section [7] deals with the details of implementing 
SDA for this problem, and in particular with the choice of the parameter 7 in the 
Cayley transform. We present several numerical experiments to illustrate the benefits 
of this approach in Section [S] and presents our conclusions in Section [9] 

2. Control and Matrix Theoretic Preliminaries. Throughout the paper 
real and complex numbers are denoted by R and C, the open left and right half-planes 
by C~ and C''", respectively. The symbol i stands for the imaginary unit, iR denotes 
the imaginary axis and by z we denote the conjugate transpose of z G C. Natural 
numbers excluding and including are denoted by N and No, respectively. The spaces 
of n X m complex matrices are denoted by C"'™, and the set of invertible and complex 
nxn matrices by G1„(C). The matrices A'^ and A* denote, respectively, the transpose 
and the conjugate transpose of A e C"'™, and A''^ = (A"!)^, A'* = (A-^)*. We 
denote by rank(A) the rank, by imA the image, by kei A the kernel, by cr(A) the 
spectrum of a matrix A. For Hermitian matrices P,Q^ C"'", we write P > Q 
{P > Q) ii P — Q IS positive (semi-)definite. 

For a rational matrix- valued function $ : C\D — > C"'™, where D C C is the finite set 



s£ -A 





si + A* 
B* 



si + A B 
Q C 
C* R 



(1.4) 
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of poles, we define the normal rank by normalrank $ = max^gcxr* rank<I>(s). 
With Ai G C"*'™' with mi,ni S Nq for i = 1, . . . we denote the block diagonal 
matrix by diag(Ai, . . . , Af^). An identity matrix of order n is denoted by /„ or simply 
by I. The zero n x m (n x n) matrix is denoted by On,m (resp. 0„) or simply by 0. 
Moreover, for fc G N wc introduce the following special matrices Jk, Mk, Nk € 
Kk,Lk G K*^-!'*^ with 



pk.k 



1 



1 01 



1 



ro 1 



1 



1 



1 



Definition 2.1. Let sE — A he a matrix pencil with E,Ae R™^". Then sE — A 
is called regular if m ~ n and normalrank(s-B — A) = n. 

A pencil sE — A is called even if E = —E* and A = A* . A pencil with E, A G 
R2"'2n ^g^iigg syniplectic if EJE'^ = AjA^ , where 



J = 



/„ 

-/„ 



Many properties of a matrix pencil can be characterized in terms of the Kronecker 
canonical form (KCF). 



Type 


Size 




Parameters 


Wl 


kj X kj 


(s - X)h^ - Nk^ 


kjGN^XeC 


W2 


kj X kj 


sNk, - Ik, 


kj e N 


W3 


i^kj 1 ^ X kj 


sKk, - Lk, 


kj e N 


W4 


kj X i^kj ' 1 ^ 


sK^ - L^ 


kj e N 



Table 2.1 

Block types in Kronecker canonical form 



Theorem 2.2. fTSf For a matrix pencil sE - A with E,Ae C"'™, there exist 
matrices Ui G G1„(C), Ur G G1,„(C), such that 

UiisE ~ A)Ur = diag(Ci(s), . . . ,Cfc(s)), (2.1) 

where each of the pencils Cj{s) is of one of the types presented in Table [KT\ 
The numbers A appearing in the blocks of type Wl are called the (generalized) eigen- 
values of sE—A. Blocks of type W2 are said to be corresponding to infinite eigenvalues. 

A special modification of the KCF for even matrix pencils, the so-called even Kro- 
necker canonical form (EKCF) is presented in[38j. 
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Type 


Size 




Parameters 


El 






Qk,.k, {\-s)Ik-Nk^ 
_{\+s)h~Nl Qu,.k, \ 


kj e N, A e C+ 


E2 


kj X kj 




kj g N, A* e K, 


E3 


kj X kj 


ej{isMk^ + Jkj) 


kj e N, 
£{-1,1} 


E4 


(2%-l)x 
(2%-l) 




sK^^+I^i^ Ok,.k, \ 


kj G N 



Block types in even Kronecker canonical form 



Theorem 2.3. For an even matrix pencil sE — A with E,A e C"'", there 
exists a matrix U G G1„(C) such that 

U*{sE - A)U = diag(Pi(5), . . .,Vu{s)), (2.2) 

where each of the pencils 'Dj{s) is of one of the types presented in Table [KM 
The numbers Sj in the blocks of type E2 and E3 are called the block signatures. 
The blocks of type El contain pairs pairs (A, —A) of generalized eigenvalues. The 
blocks of type E2 and E3 respectively correspond to the purely imaginary and infinite 
eigenvalues. Blocks of type E4 consist of a combination of blocks that arc equivalent 
to those of type W3 and W4. 

Definition 2.4. A subspace V C is called (right) deflating subspace for the 
pencil sE — A with E,A^ {£M,n ij jor a matrix V G C^'*^ with full column rank and 
\TD.V — V, there exists an I < k and matrices W G C^^"', E,Ag with 

{sE - A)V = W{sE - A), (2.3) 



Definition 2.5. An eigenvalue X of a matrix pencil is called c-stable, c-critical 
or c- unstable respectively «/Re(A) is smaller than, equal to, or greater than 0. A right 
deflating subspace is called c-stable (resp. c- unstable^ if it contains only c-stable (resp. 
c-unstable) eigenvalues, and c-semi-stable (resp. c-semi-unstable^ if it contains only 
c-stable or c-critical (resp. c-unstable or c-critical) eigenvalues. The same definitions 
hold replacing the prefix c- with d- if we replace the expression Re(A) with |A| — 1. 

Definition 2.6. Let M G C*^-*^ be given. A subspace V dC^ is called Al-neutral 
if x*My = for all x, y eV. 

Definition 2.7. Let a pair iA,B) G C"-" x C"'™ be given. Then 
(i) {A,B) is called controllable if xank[sL — A , B] = n for all s G C; 
(ii) {A,B) is called stabilizable i/ rank[ s/— A , B] ^ n for all s G C+. 

Definition 2.8. Given 7 G M, 7 7^ 0, the Caylcy transform of a regular pencil 
s£ — A is the pencil 
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The Cayley transform preserves left and right eigenvectors and Jordan/Kronecker 
chains, while transforms the associated eigenvalues according to the map C : A i— >• j^, 
A S C U oo. In particular, Kronecker blocks of size k for A are mapped to Kronccker 
blocks of size k for C(A). 

3. Solvability of Lur'e equations. In this part we collect theoretical results 
being equivalent for the solvability of Lur'e equations. For convenience, we will call 
a Hermitian matrix X a solution of the Lur'e equations if p.ip is fulfilled for some 
K e CP'", L e CP'™. 

We now introduce some further concepts which are used to characterize solvability 
of the Lur'e equations. 

Definition 3.1. For Lur'e equations the spectral density function is de- 

fined as 



'{iLuI ~A)-^B 



- * 




c 




{iiol 




c* 


R 







(3.1) 



In several works, the spectral density function is also known as Popov function. 
Definition 3.2. For Lur'e equations U.l]) . the associated linear matrix inequal- 
ity (LMI) is defined as 

'A*Y + YA + Q YB + C 
B*Y + C* R 

The solution set of the LMI is defined as 

Slmi = {Y e C"'" : Y is Hermitian and ( TO)] holds true}. (3.3) 

The LMI 113.2]} is called feasible if Slmi 0- It can be readily verified that Y E Slmi 
solves the Lur'e equations if it minimizes the rank of p.2[) . We now collect some known 
equivalent solvability criteria Lur'e equations. In the following we require that the 
pair (A, B) is stabilizablc. Note that this assumption can be further weakened by 
reducing it to sign- controllability |35j . This is not considered here in more detail. 

Theorem 3.3. Let the Lur'e equations with associated even matrix pencil 

s£ ~ A as in |_?./^[ ) and spectral density function $ as in iS. 1\) be given. Assume that 
at least one of the claims 

(i) the pair [A, B) is stabilizablc and the pencil s£ — A as in ^1.4^ is regular; 

(a) the pair ( A, B) is controllable; 
holds true. Then the following statements are equivalent: 

1. There exists a solution {X,K,L) of the Lur'e equations. 

2. The LMI is feasible 

3. For all oj eM. with iuj ^ a{A) holds $(ia;) > 0; 

4. In the EKCF of s£ — A, all blocks of type E2 have positive signature and even 
size, and all blocks of type E3 have negative sign and odd size. 

5. In the EKCF of s£ — A, all blocks of type E2 have even size, and all blocks 
of type E3 have negative sign and odd size. 

In particular, solutions of the Lur'e equations fulfill {X,K,L) G C"'" x C"'P x C™'^ 
with p = normalrank $ . 



> 0. 



(3.2) 



It is shown in [35] that m — normalrank $ is to the number of blocks of type E4 
in an EKCF of s£ — A. In particular, the pencil s£ — A is regular if and only if $ has 
full normal rank. 
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x+ 





In 











Now we place particular emphasis on the so-called maximal solution. 

Theorem 3.4. Let the Lur'e equations lll.l]) be given with stabilizable pair {A, B). 
Assume that Slmi o,s defined in Ii3.3fl is non-empty. Then there exists a solution 
of the Lur'e equations that is maximal in the sense that for all Y G Slmi holds 

Y < X+. 

The following result [3S] states that the maximal solution can be constructed via the 
c-stable deflating subspace of the associated even matrix pencil s£ — A. 

Theorem 3.5. Let the Lur'e equations U.l]) be given with stabilizable pair {A, B). 
Assume that Slmi cis defined in US. 3\) is non-empty. Then 



(3.4) 



spans the unique {n + m)- dimensional semi-c-stable £ ^neutral subspace of the pencil 
(|1.4p . The above theorem states that the maximal solution can be expressed in 
terms of a special deflating subspace of s£ — A. By means of the EKCF, this space 
can be constructed from the matrix U G (J32n+m,2n+m ^jj-ji^ging the pencil s£ — A into 
even Kronecker form (|2.2p . Considering the partitioning U = [Ui , . . . , Uk] according 
to the block structure of the EKCF, a matrix V £ £^'2n+m,n+m spanning the desired 
deflating subspace can be constructed by 

V^[Vi ... Vk] for Vj = UjZ.j, (3.5) 

where 

[/fe^,Ofc^. if is of type El, 

2. ^ } ' if T^j is of type E2, 

^ ' [I(k,+i)/2 , 0(fc^_i)/2]'^, if 'Dj is of type E3, 
hj , Ofc^.+i f, if Vj is of type E4. 

In particular, the desired subspace contains all the vectors belonging to the Kronecker 
chains relative to c-stable eigenvalues, no vectors from the Kronecker chains relative 
to c- unstable eigenvalues, the first kj/2 vectors from the chains relative to c-critical 
eigenvalues, and the first (kj + l)/2 from the chains relative to eigenvalues at infinity. 

4. Outline of SDA. The structured doubling algorithm (SDA) is a matrix it- 
eration which computes two special deflating subspaces of a matrix pencil. It was 
introduced by Anderson [4] as an algorithm for the solution of a discrete-time alge- 
braic Riccati equation, and later adapted to many other equations and explained in 
terms of matrix pencils in several papers by Wen- Wei Lin and others [TTJ [521 US [TO] . 
It is strongly related to the sign function method and to the disc method for matrix 
pencils [Sl[7]. 

Theorem 4.1 ([S]). LetA-s£ be a regular matrix pencil with A, £ G ^n-{^mm-\-m ^ 
and let A^t — s£^ be a regular pencil of the same size with A*£ — £tA. Then 

1. the pencil ArA — sf*f is regular and has the same right deflating .subspaces 
as A — s£ 

2. its eigenvalues are the square of the eigenvalues of the original pencil. 
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This result is far from surprising in the case in which £ is invertiblc: in this case, 
the eigenvalues and right deflating subspaces oi A — s£ correspond to the eigenvalues 
and right invariant subspaces oi£~'^A. and it is simple to check that the conditions im- 
posed on A, imply {£^£)-'^{A^A) = {£~'^Af. Thus the map {A,£) ^ {A^A,£^£) 
is a way to extend the concept of squaring to matrix pencils. 

A pencil A — s£ with A,£€ '^'^+m,n+m jg ^^^^^ -^^^ standard symplectic-like 

form I (SSF-I) if it can be written as 



A^ 



E 
-H Im 



£ ^ 



In 




-G 
F 



(4.1) 



where the block sizes are such that E G K^^^ and F g R^^^*^. Note that a pencil 
in SSF-I is always regular. When a pencil A — s£ is in SSF-I, a choice of — s£^, 
satisfying the requirements of Theorem 14.11 is 



E{In - GH)-^ 
-F{Im-HG)-^H 





In 



In 




-E{In -GH)-^G 
F{Im-HG)-^ 



yielding a new pencil A^,A — s£^£ which is still in SSF-I: 



A = 



E{In -GH)-^E 
-{H + F{Im-HG)-^HE) In 



£ = 



In 




-{G + E{In -GH)-^GF) 
F{Im - HG)~^F 



The only hypothesis needed here is that / — GH and / — HG are nonsingular. In fact, 
by the Sherman-Morrison formula, they are either both singular or both nonsingular. 
We may design Algorithm [T] by repeating this transformation. Each step of the 



input : Eq, Fq, Go, Hq defining a pencil in SSF-I 

output: HoojGoo so that the subspaces in (|4.2p are respectively the canonical 
semi-d-stable and semi- d- unstable deflating subspaces of the given 
pencil 

k< — 0; 

while a suitable stopping criterion is not satisfied do 

E* < — EkilN — GkHk)^^; 

F* < — Fkihi — IIi;Gk)~^; 

Gk+i < — Gk + E^GkFk] 

Hk+i < — Hk + F^HkEk] 

Ek+i < — E^Ek] 

Fk+i < — F^Fk; 

k i — fc + 1; 
end 

Hod —Hk', 
Goo —Gk', 

Algorithm 1: SDA-I 



algorithm costs ^{M"^ + N"^) + 6MN{M + N) floating point operations. This reduces 
to in the case M = N. 

The following result is proved in for the symplectic case and in [T^ for several 
specific matrix equations, but its proof works without changes for our slightly more 
general case. 

We shall call a pencil weakly d-split if there exists an r such that: 
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• the lengths of the Kronccker chains relative to d-stable eigenvalues sum up 
to TV - r; 

• the lengths of the Kronecker chains relative to d-unstable eigenvalues sum up 
to M - r; 

• the lengths kj of the Kronecker chains relative to d-critical eigenvalues (which 
must sum up to 2r if the two previous properties hold) are all even. 

In this case, we define the canonical d-semi-stable (resp. d- semi- unstable) subspace as 
the invariant subspace spanned by all the Kronecker chains relative to d-stable (resp. 
d-unstable) eigenvalues, plus the first kj/2 vectors from each critical chain. 

Theorem 4.2. Let the pencil (j4.ip be weakly d-split, and suppose that there are 
two matrices in the form 



'In' 




Goo 









(4.2) 



spanning respectively the canonical d-stable and d-unstable deflating subspace. Then 
for Algorithm{l\ it holds that 

. \\E,\\=0{2-'^), 

. ||F,|H0(2-'=), 

. \\Hoo-Hk\\=0{2-''), 

• IIGoo-Gfcll =0(2-^). 
Notice that, when N = M, a pencil in SSF-I is symplectic if and only if E* — F, 
G = G* and H = H* . In this case, all the pencils generated by the successive steps of 
SDA are symplectic, i.e., at each step k we have = Fk, Gk = G^, Hk = H^. The 
implementation can be slightly simplified, since there is no need to compute Ek+i 
and -Ffc+i separately, nor to invert both In — GkHk and Im — HkGk, as the second 
matrix of both pairs is the transposed of the first. 

5. A method to compute the SSF-I of a pencil. One can transform a 
regular pencil into SSF-I easily using the following result. 

Theorem 5.1. Let A - s£ be a matrix pencil with A,£ e ^n+m.n+m ^ 
partition both matrices as 

A = [^1 A2] £ = [£i £2] 

with Ai,£i e and A2,£2 e R^+^-^'^^ A SSF-I pencil having the same 

eigenvalues and right deflating subspaces of the original pencil exists if and only if 

[£i A2] (5.1) 

is nonsingular; in this case, it holds 

= [£i A2]~'[Ai £2]. (5.2) 



E -G 
-H F 



Proof. We are looking for a matrix Q such that 

g [^1 A2] - sQ 1 £2] 

By taking only some of the blocks from the above equation, we get 



" E 





— s 


I 


-G 


-H 


I 





F 



Q£i = 



, QA2 = 




/ 
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I.e. 



I 
/ 



thus Q must be the inverse of the matrix in (|5.1I) . 

On the other hand, taking the other two blocks we get 



E 

-H 



, Q£2 = 



-G 

F 



which promptly yields (j5.2p . □ This formula is strictly related to the principal pivot 
transform (PPT) gD]. 

Wc point out an interesting application of Theorem 15. 1[ which is not related to 
the rest of the paper. SDA is often used in the solution of nonsymmetric algebraic 
Riccati equations (NARE) [22], where it is applied to the Cayley transform of the 
matrix 



n = 



D 
B 



with B, C, D blocks of suitable size associated with the coefficients of the problem. 
In pencil form, its Cayley transform given by (T^ — 7/) — s{TL + jI), thus (|5.2p becomes 



" E 


-G 




-H 


F 





D + 7/ -G ' 
B -A-jI 



D - 7/ -G ' 
B -A + jI 



(5.3) 



This formula is more compact to write and more computationally effective than the 
one suggested by Guo et al. [13]. In fact, their expressions for the starting blocks 
involve computing the inverses of two N x N and M x M matrices, which are indeed 
the (1,1) and (2,2) blocks of the matrix to be inverted in (|5.3p and their Schur 
complements. Clearly, two these inversions are redundant in (|5.3p . which requires 
more or less half of the computational cost with respect to the original formulas in 
[23]. 



6. A reduced Lur'e pencil. Let s£ — A be the pencil in (|1.4|) . By Theorem l5.1 
the SSF-I form of its Cayley transform (assuming M = n + m, N = n) is given by 



E 



Let now 



-G 
F 



A~-fI B 
A* --fl Q G 
B* G* R 



A 



A -7/ 
A* - 7/ Q 



A + "fI B 
A* +jl Q G 
B* G* R 



,S 



/„ 



and assume that both A and its Schur complement 



R - [B* G*] A'^ 



B 
G 



= $(7) 
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(here $ is the same function as in p.ip , as one can verify by expanding both expres- 
sions) are nonsingular. In this case we can perform the inversion with the help of a 
block LDU factorization. Tedious computations lead to a matrix of the form 



A 

B Im. 



with 



/ + 27^-^5' + 27 A" 1 



$(7)-! [B* C*] A-^S, S 



In the blocks used in SSF-I, this means that 



/ 
/ 



(6.2) 



F 



F 



G 



H 



where the blocks F, G, H have size n x n. It follows that a special right deflating 
subspace of this pencil is 



02nxm 



whose only eigenvalue is 1 with algebraic and geometric multiplicity m, while the 
deflating subspaces relative to the other eigenvalues are in the form 

V' 
* 

where V has 2n rows and is a deflating subspace of the reduced pencil 





-G 




■ E 


0" 





F 




-H 


In. 



(6.3) 



Using (|6.2p and the fact that A and $(7) are symmetric, one sees that AS is sym- 
metric, too. This means that E* = F and G = G*, H = H*, that is, the pencil 
is symplectic. 

The pencil (|6.3p is given by P*{s£ — A)P, where P is the projection on 



span 



(ker£:)- 



With this characterization, it is easy to derive the KCF of (|6.3|) from that of the 
Cay ley transform of (|1.4p . We see that kerf is the space spanned by the first column 
of each W2 block (as a corollary, we see that there are exactly m = dim ker E such 
blocks). These blocks are transformed into blocks Wl with A = 1 by the Cayley 
transform. Thus projecting on their orthogonal complement corresponds to dropping 
the first row and column from each of the Wl blocks relative to A = 1. In particular, 
it follows that if the criteria in Theorem 13.31 hold, then every Kronecker block of (|6.3p 
relative to an eigenvalue A on the unit circle has even size. Therefore, the reduced 
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pencil (|6.3p is weakly d-split. By considering which vectors are needed from each 
Kronecker chain to form the subspace in p.5|) . we get therefore the following result. 

Theorem 6.1. LetV be the unique {n + m)- dimensional c-semi-stable £-neutral 
deflating subspace of \1A\ . Then, there is a matrix V2 G C"'™ such that 



V = span 



V2 







where Vi spans the canonical d- semi-unstable subspace oj the pencil (j6.3p . Moreover, 
i/span(Vi) admits a basis in the form 

In 



_l_ is the maximal solution of the Lur'e equation (jl.ip . 



then X 

the canonical weakly stabilizing solution of the DARE 

-1^, + 



X = EX{I ~ HX) 



G. 



In other words, X is 



(6.4) 



7. Implementation of SDA. Based on the results of the previous sections, 
we can use the SDA-I algorithm to compute the solution to a Lur'e equation. The 
resulting algorithm is reported as Algorithmic] 

As we saw in Section 21 the symplecticity of the pencil is preserved during the 
SDA iterations, and helps reducing the computational cost of the iteration. Moreover, 
this has the additional feature of preserving the eigenvalue symmetry of the original 
pencil along the iteration. 

The explicit computation of (a possible choice of) K and L are typically not 
needed in the applications of the Lur'e equations. If they are needed, they can be 
computed using the fact that 



A*X + X*A + Q XB + C 
B*X* + C*R 





'K*' 




L* 



[K L] 



(7.1) 



is a full rank decomposition. 

The accuracy of the computed solution depends also on an appropriate choice of 
7. Clearly, two goals have to be considered in the choice: 

1. the matrix to invert in (|6.ip should be well-conditioned; 

2. the condition number of the resulting problem, i.e, the separation between the 
stable and unstable subspace of the Cayley-transformed pencil (|6.3p . should 
not be too small. 

While the impact of the first factor is easy to measure, the second one poses a more 
significant problem, since there are no a priori estimates for the conditioning of a 
subspace separation problem. Nevertheless, we may try to understand how the choice 
of the parameter 7 of the Cayley transform affects the conditioning. Roughly speak- 
ing, the conditioning of the invariant subspace depends on the distance between its 
eigenvalues and those of the complementary subspace ^U\. The eigenvalues of the 
transformed pencil are given by = 1 — ^^p^, thus they tend to cluster around 1 
for small values of 7, which is undesirable. The closest to 1 is the one for which A -I- 7 
has the largest modulus; we may take as a crude approximation p{A) + 7, which can 
be further approximated loosely with \\A\\^ -I- 7. Therefore, as a heuristic to choose a 
reasonable value of 7, we may look for a small value of 

\A\ 



7(7) = max I condest ( [£i A2] ) 



II 1+7 
27 
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input : A, B, C, Q, R defining a Lur'e equation 

output: The weakly stabilizing solution X (and optionally K and L) 

Choose a suitable 7 > 0; 

Compute 



T 



A-^I B 
A* -jl Q C 
B* C* R 



A + 7/ 
A* + 7/ Q 
B* C* 



Partition 



T = 



E -G 
-H E* 



Use SDA-I on E', F = E* , G,H to compute Goo, Hoo, and set X = Goo; 
if K and L are needed then 

Compute S and V* corresponding to the first m singular vectors of the 
SVD of JLI]); 
Set [A' L] i — S;i/2y*; 
end 

Algorithm 2: A SDA algorithm for the maximal solution of a Lur'e equation 



where condest(-) is the condition number estimate given by Matlab. Following the 
strategy of [11] , we chose to make five steps of the golden section search method [30] 
on 7(7) in order to get a reasonably good value of the objective function without 
devoting too much time to this ancillary computation. 

However, we point out that in |11| . a different /{"j) is used, which apparently only 
takes into account the first of our two goals. The choice of the function is based on an 
error analysis of their formulas for the starting blocks of SDA. Due to Theorem 15. 2[ 
this part of the error analysis can be simplified to checking condest ( [£i A2] ) . 

8. Numerical experiments. We have implemented Algorithmic] (SDA-L) using 
Matlab, and tested it on the following test problems. 

PI a Lur'e equation with a random stable matrix A G R"'", a random C = B, Q = 
and R the m x m matrix with all the entries equal to 1, with rank(i?) = 1. 
Namely, B was generated with the Matlab command B=rand(n,m), while to 
ensure a stable A a we used the following more complex sequence of com- 
mands: V=:randn(n); W=randn(n); A=-V*V'-W+W'; 

P2 a set of problems motivated from real-world examples, taken with some modifica- 
tions from the benchmark set carex [8] . Namely, we took Examples 3 to 6 (the 
real-world applicative problems) of this paper, which are a set of real-world 
problems varying in size and numerical characteristics, and changed the value 
of R to get a singular problem. In the original versions of all examples, R 
is the identity matrix of appropriate size; we simply replaced its (1,1) entry 
with 0, in order to get a singular problem. 

P3 a highly ill-conditioned high- index problem with m = 1, A = In + Nn, B ~ Cn 
(the last vector of the canonical basis for M"), C = —B, R = 0, Q = 
— tridiag(l, 2,1). Such a problem corresponds to a Kronecker chain of length 
2n+l associated to an infinite eigenvalue, and its canonical semi-stable solu- 
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Fig. 8.1. Relative residual for PI 
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Fig. 8.2. Relative residual for P2 
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tion is X — I. Notice that the conditioning of the invariant subspace problem 
in this case is ei/(2"+i)^ for an unstructured perturbation of the input data 
of the order of the machine precision e [IHl section 16.5]. 
The results of SDA-L are compared to those of a regularization method as the 
one described in (|1.3|) . for different values of the regularization parameter e. After the 
regularization, the equations are solved using Algorithm [T] after a Cayley transform 
with the same parameter 7 (R+S), or with the matrix sign method with determinant 
scaling [32j [24] (R+N) . We point out that the control toolbox of Matlab contains a 
command gcare that solves a so-called generalized continuous-time algebraic Riccati 
equation based on a pencil in a form apparently equivalent to that in (|1.4p . In 
fact, this command is not designed to deal with a singular R, nor with eigenvalues 
numerically on the imaginary axis. Therefore, when applied to nearly all the following 
experiments, this command fails reporting the presence of eigenvalues too close to the 
imaginary axis. 

For the problem P3, where an analytical solution X ~ I is known, we reported 
the values of the forward error 

X -X 

For PI and P2, for which no analytical solution is available, we computed instead the 
relative residual of the Lur'c equations in matrix form 

lA'X + XA + Q XB + C 
[ B*X* + C* R 

\A'X + XA + Q XB + Cl 

B*X*+C* R 

L -1 r 

(see ()7.ip ). A star ★ in the data denotes convergence failure. 

We see that in all the experiments our solution method obtains a better result 
than the ones based on regularization. 

9. Conclusion and open issues. In this work we have introduced a new nu- 
merical method for the solution of Lur'e matrix equations. Unlike previous methods 
based on regularization, this approach allows one to solve the original equation with- 
out introducing any artificious perturbation. 





K* 




L* 
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Fig. 8.3. Forward error for PS 
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The first step of this approach is applying a Cayley transform to convert the 
problem to an equivalent discrete-time pencil. In this new form, the infinite eigenval- 
ues can be easily deflated, reducing the problem to a discrete-time algebraic Riccati 
equation with eigenvalues on the unit circle. For the solution of this latter equation, 
the structured-preserving doubling algorithm was chosen, due to its good behaviour 
in presence of eigenvalues on the unit circle, as proved by the convergence results in 
[25] , Direct methods, such as the symplectic eigensolvers presented in [17], could also 
be used for the solution of the deflated DARE. 

The numerical experiments confirm the effectiveness of our new approach for 
regular matrix pencils. It is not clear whether the same method can be adapted to 
work in cases in which the pencil (|1.4p is singular, which may indeed happen in the 
contest of Lur'e equations. Another issue is finding a method to exploit the low-rank 
structure of Q (when present). These further developments are currently under our 
investigation. 
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